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^ We introduce a fast Fourier transform on regular J-dimensional lattices. 

We investigate properties of congruence class representants, i.e. their order- 
ing, to classify directions and derive a Cooley-Tukey-Algorithm. Despite 
the fast Fourier techniques itself, there is also the advantage of this trans- 
^ form to be parallelized efficiently, yielding faster versions than the one- 

^ dimensional Fourier transform. These properties of the lattice can further 

be used to perform a fast multivariate wavelet decomposition, where the 
wavelets are given as trigonometric polynomials. Furthermore the preferred 
f^ directions of the decomposition itself can be characterised. 

fS| Keywords, wavelets, lattices, multivariate fast Fourier transform, periodic multireso- 

>^ lution analysis, Dirichlet wavelets, shift invariant space 

^ 1 Introduction 

^ Recently, a framework for multivariate periodic wavelet analysis [15] was developed, 

^H using shift invariant spaces of translates defined by a regular integral matrix M. We 

>. further investigate this periodic multiscale analysis and develop fast algorithms for the 

'k> decomposition of a periodic multivariate function. This generalises the one-dimensional 

Vh case of periodic wavelets as described e.g. in [18, 20, 21]. 

There are many approaches towards decompositions of multivariate functions in terms 
of wavelets that have been investigated in the last two decades, e.g. curvelets [5, 6], 
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ridgelets [4], contourlets [9] or shearlets [8]. The periodic wavelet transform is defined 
on a finite set of translates and hence — in contrast to the forementioned wavelets — leads 
to finite sums by construction. These translates are defined by the pattern of a matrix M 
and lead to a multivariate Fourier transform onto a frequency domain with the shape of 
a parallelogram. The corresponding multivariate Fourier matrix is a Kronecker product 
of one-dimensional Fourier matrices, as mentioned by [1]. 

The main challenge for an implementation is the ordering of the elements. While the 
elements on the real line are naturally ordered, there are many different ways to order the 
points of a lattice. We combine the idea of using the Smith normal form, as described 
by Mersereau et al. in [16, 17] with the theoretical results on abelian groups developed 
by Auslander et al. in [1, 2]. While the former authors use arbitrary decompositions 
of the regular matrix M to derive a multivariate Cooley-Tukey-Algorithm, we use the 
approach from the latter authors and split the pattern into subgroups. Furthermore, our 
investigation also improves the algorithm mentioned in [16] by using an order in terms 
of specific bases. 

Despite the mathematical motivation as a generalisation of the rectangular multivari- 
ate Fourier transform, the lattice based Fourier transform naturally arises in crystallog- 
raphy. Based on X-ray defractions, the shape of the so called unit cell and its internal 
structure of a crystal is examined, e.g. for proteins [10]. The reciprocal lattice used 
for the measurements of X-Ray defractions corresponds to the generating group from 
Section 2. The reciprocal lattice is derived by looking at biorthognal bases [12, Chapter 
2], which are also used here in Section 3. 

The paper starts with some preliminary investigations and notations in Section 2. In 
Section 3, we define and characterise bases for the patterns and use these to address all 
elements. This also leads to the notion of the dimension of the pattern and describes an 
ordering with respect to specific bases. Using these orderings, there is no rearrangement 
necessary for the Fourier transform, neither in time nor in frequency. This is used in 
Section 4 to develop a fast algorithm for the Fourier transform on the pattern that can 
additionally be parallelized. In Section 5, we describe properties of bases, especially 
together with the bases of the dual group — the so called generating group of a matrix 
M. These are applied to the multivariate periodic wavelet transform developed in [15] to 
obtain a fast algorithm for the decomposition. Furthermore, we are able to use the basis 
of the pattern and generating group to generalise the one-dimensional scaling property 
given in [21], including the characterisation of directions and their transformation from 
one scaling space to the next. 

Finally in Section 6, we discuss an example for both the Fourier transform and the 
wavelet transform. The first example is concentrated on computational costs of the 
multivariate Fourier transform and possible parallelizations. The wavelet transform is 
presented by decomposing two different box splines. 



2 Preliminary investigations 
2.1 Function spaces 

The space of functions under consideration is the Hilbert space L?{T'^) of all square 
integrable functions on the torus T'' = [0, 2nY with the inner product 



{f,g) = j^l^/{x)g{x)dx, for/,gGL2(T^). (1) 



Every function / G L?{T'^) can be decomposed with respect to the monomials e'"* ° E 
L^{T^), k G Z^, which yields the Fourier series representation 

/(x) = £ Ck(/)e^'^^ where Ck(/) = (/,e*^°), k G Z^ (2) 

We denote by c(/) = {c\i{f))^^^d G 1^{Z^) the generalised sequences which form a 
Hilbert space with the inner product 

(c,d)= £ CkJ,^, c,dG/2(Z^), 
keZ'' 

where the Parseval equation reads as 



{f,8) = W)Ms))= Ick(/)ck(^). 

keZ'' 

We define for any y G M'^ the translation operator T{y)f = /(o — 2;ry), / G L^{T^). A 
straight forward computation leads to Ck(r(y)/) = e^^^^^ ^Ck(/). We can restrict y to 
any shifted unit cube, e.g. [—j, jY, because q^"^^^^ ^ = 1, k,z G Z'^. 

2.2 The pattern and the generating group 

For any regular matrix M G Z'^^'^, we define the congruence relation for h, k G Z"^ with 
respect to M by 

h = k mod M ^ 3z G Z'^ : k = h + Mz. 

We define the lattice 

A(M) := M^iZ'^ = {ye W' : My G Z^}, 

note that it is 1 -periodic and define the pattern ^(M) as any complete set of congruence 
class representants, e.g. A(M) fl [0, lY or A(M) H [— ^ , ^) , which both contain exactly 



one element of each congruence class with respect to mod I on A(M). We denote by 
[x]m the congruence class of x G A(M) and define by 

x|^(M):= WMn=^(M), xGA(M) 

the mapping from any point x G A(M) of the lattice onto its congruence class represen- 
tant belonging to ^(M). Then (^(M), +|_^(m)) is an abelian group. 

Since M G Z^^'^ is regular, it can be seen as a bijective map. Hence all the definitions 
and properties above also hold for the generating group ^(M) := M,^(M) equipped 
with +|^(M)' where Mo : ^(M) — )■ ^(M) performs a group isomorphism between 
(^(M), +|.^(M)) and (=^(M), +|^(m))- 

The Smith normal form is defined by the decomposition 



M 



vdxd „,u„„„ "c a:^^ ( r. \" 



QER, Q,E,RgZ^^'', where E = diag(e/)^j , (3) 

with |detR| = |detQ| = 1 and the elementary divisors £j G N — which exist due to 
the Theorem on elementary divisors, see e.g. [14, Chapter 10] — fulfill £j\£jj^i, j = 
1, . . . ,(i — 1. Besides the already mentioned isomorphism this also implies that 

^(M) = ^(E) = ^(E) = ^(M) = %^ ® ^£2 ® ■ ■ ■ ® ^firf' (4) 

where "^^ = {0, 1, . . . ,z — l},z G Z denotes the cyclic group Z/zZ. Hence |^(M) | = 
I ^ (M) I = Ci ■ . . . ■ Cj = I detM| =: m. We further denote by Jm := #{^j > 1 } the number 
of cycles greater than 1 . 

For any decomposition of a regular matrix M = JN, N, J G Z"^^"^ it holds [15], that 
there exists a unique decomposition for y G ^(M) 

y=(x + N-iz)|^(j^), xG^(N),zG^(J) (5) 

which can also be applied to ^(M^) yielding for h G ^(M^) the unique decomposition 

h= {g + N^k)l^^^ry kG^(J^), gG^(N^). 

The Fourier transform on the pattern ^(M) is defined [7] by 
^(M) = — fe-2^ih^M-lg^ _ _±_ ('e-2mh^y\ 



(6) 

where h G '^(M'^) indicates the rows, y G ^(M) indicates the columns and the equality 
holds if the ordering of the elements in ^(M) and ^(M) are identical with respect to 
the bijection Mo. The discrete Fourier transform on ^(M) is defined for a vector 
a = («y)y£^(M) ^ C'" arranged in the same ordering as the columns in (6) by 

a = (^h)heff (M^) = ^(M)a, 
where the vector a is arranged as the columns of ^(M) in (6). 
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Figure 1. The pattern ^(M) for M given in (7) is illustrated on the left, where 
the set of congruence classes is chosen from [— ^, ^)^. The emphazised point 

yi ~ (I' J2) ^^ ^^^ canoncial basis vector of the pattern from Section 3. Its 
corresponding generating group ^(M^) (right) illustrates the set of frequen- 
cies obtained when performing a Fourier transform on the points of ^(M) 
using the matrix from Eq. (6). 



Example 1. We a look at the sheared and rotated 2-dimensional lattice generated by 
the matrix 



M 



4 -3 
4 5 



1 -iW4 r 
1 1 Mo A)- 



(7) 



One way of choosing the pattern ^{M.) is illustrated in Fig. 1 on the left. The congru- 
ence class representants are chosen from [—2,2)^- They can also be seen as sampling 
points on the 1 -periodic torus. The lattice A(M) consists of all integer shifts of these 
points. The corresponding generating group ^(M-^), used to construct the Fourier ma- 
trix (6), is given in Fig. 1 on the right. They are obtained by taking all integer valued 



vectors from M 



1 ]_\d 

'2' 2' 



The Smith normal form is given as 



M 



3 1 

5 2 



1 

32 



12 



QER, 



hence both groups are isomorphic to ^32. 



3 Bases for the pattern and the generating group 

In contrast to the one-dimensional pattern, i.e. the set {0, \/N, ...,(//— 1)/A^}, the or- 
dering of the elements in ^(M) and ^{M^) is not given by any "natural" topology. 
Due to the isomorphisms in (4) this section deals with finding an ordering, that is sim- 
ilar to the tensor product case by introducing a basis for ^(M). Multiplying all the 
formulae in the following section with M leads to the same properties for the generating 
group ^(M), hence starting with ^(M^) also a basis for '^(M^) can be constructed. 
For a fixed regular matrix M we define the set of vectors 

y;:=R"^- ^d-du+j^ 7 = 1,...,Jm, (8) 

where R denotes the left basis transform in the Smith normal form (3) and ej denotes the 
7-th unit vector. These vectors are linear independent, because R has full rank. Using R 
as a change of basis, we see from [15, Lemma 2.4] that A(M) = A(ER) and 

A(M) = A(ER) = {y | ERy G Z''} = A(E) = {x | Ex G Z^}, 

where / : A(ER) — )■ A(E),/(y) = Ry is an isomorphism between A(M) and A(E). 

The scaled unit vectors e^}^ ^ j^d-dM+ /' J = 1 ; • • • ; ^m, form a basis for ^ (E) , i.e. 
there is a unique representation for each x G ^(E) of the form 



L ^j ed-,^+j^d-dM- 



, where < A^ < Cd-du+j- (9) 



The other unit vectors are not scaled due to j < d — Jm "^ ^j = 1' ^^^ hence van- 
ish in such a summation with respect to the congruence classes, i.e. modi. By us- 
ing the inverse of the isomorphism /^^ : A(E) — )■ A(ER),/^^(y) = R^^y, the vectors 

{yi. • • • .yrfM^ f°™^ ^ ^^^^^ °f <^(M). 

Remark 1. The lattice generated by {yi, . . . ,yrf^} is called rank-dM-l^ttice. Further 
the ordering of the basis elements is crucial, because each vector y j has to span a cycle 
of length £d-dM+j (with respect to -Ig^iM)) '^ this notation. The value Jm is also called 
the dimension of the pattern ^(M). 

Using the basis we obtain an ordering of the elements of ^(M) by using the lexico- 
graphical ordering of 



^(M) = 
on the set Em = {0, 1 




\ £d-dyi + l-^,--,£d-^ 




^m) (A,,...,A,^)=0 




erf-rfM+i}x---x{0,l,... 


,£d} of indices 



Example 2. For the matrix M = ( ^ ^^ ) from Ex. 1 we see, that Jm — 1 (^nd the only 

basis vector is yi = R^^e2 = 32(12 l) =(§732) > which is emphasised in Fig. 1 
(left). Hence by writing ^ (M) = { fcy j | j^/j^-j , fc = 0, 1 , . . . , 3 1 } we even get an ordinary 
ordering of the elements. 

Applying these ideas on the generating group, leads to a basis of ^(M^) denoted by 

h; := R^ed-dM+j^ j = 1, . . . , 6?m. (10) 

This leads to the biorthogonality of these two bases, which is shown in the following 

Lemma 1. Let Me Z'^^'^ be regular Then the bases of^iM) and'^iM^) given in (8) 
and (10) are biorthogonal, more precisely 



1 



ifi = i, 



V/,j €{!,..., Jm}: (h;,y,) = <^'''-''M+'- (U) 

[0 else. 

Proof. For arbitrary /, j G { 1 , . . . , Jm} it holds 

(hj,y/) = (R erf_rf^j+;,R" ed_rf^j+/) = ^d-dM+fd-dm+i- 

'^d—dM+i ^d—dM+i 

n 



4 The fast Fourier transform on ^(M) 

This section is devoted to derive a fast Fourier transform on ^ (M) using its charac- 
terisation obtained in the previous section. The idea generalises an algorithm described 
in [16, 17]. It does not depend on finding a prime factor decomposition of M as men- 
tioned in [17, Eq. (35)]. It uses the Smith normal form similar to [16], but our approach 
is able to omit the rearrangement steps. 

A more general approach for abelian groups and their characters can be found in [1]. 
In contrast, our approach uses properties of the pattern and generating group to deduce 
a fast algorithm working on arrays. 

4.1 Properties of the multivariate Fourier transform 

Based on [15, Lemma 2.1], there exist permutation matrices Ph,Py on ^(M-^) and 
^(M) respectively, such that 

^(M)=Ph(#-e,®#-£2®---®=^ejPy, (12) 



where 



^. 



v^ 






eeW 



denote the elementary divisors of the Smith normal form of M. Using the bases con- 
structed in Section 3, this section develops properties which simplify (12) to 



^(M)=^ 



^ 



)^ 



)^ 



£d-dM + l ^ ^ £d-dM+2 ^ ^ -^ Erf' 



(13) 



where the first factors JF^^ , . . . , ^e,;_^ of the Kronecker product in (12) can be omitted 
due to ^1 = e'^ = 1 . The following theorem characterises conditions for the ordering 
which are necessary for ^(M) to fulfill (13). We start with the orderings using the bases 
from (8) and (10) of ^(M) and '^(M^) from Section 3 and look at their lexicographical 
orderings, i.e. 



^{M) 




and ^{M' 



where A = (Ai, . . . , A^^) and ii = {iiu.. . ,jUrf^j). 




' (14) 



f^ 



(M^)/ M^Em 



Theorem 1 (A basis for the Kronecker product). 

The orderings from (14) used for the construction of JF (M.) fulfill (13). 

Proof. Given any j E {d — Jm + 1, ■ • • ,^}, in particular £j > 1, orderings in (14) and 
Lemma 1 are used to prove validity of (13) by induction over Jm- For Jm = 1 the 
matrix ^(M) is identical to the one-dimensional case, hence Lemma 1 implies (13) . 

For Jm > 1 let Ai , . . . , A^,^-! and /ii , . . . , jU^m-i be given, each fulfilling < Aj, /i^ < 
£d-dM+j^ J = 1 , . . . , Jm, and denote 



1 



rfM-l 



Then the submatrix of ^(M) induced by the elements 



x/ 



and 



Zfc 






dM 
7=1 



^ = ^du = O5 • • • , £rf - 1 



^(M) 



^ = MrfM = 0, . . . , Crf - 1 



^(M?') 



can be written as 



,-2mzi x/ 



kj=0 



'M 



'M 



-2m\^^^iJhJ [L^jyj 



V 



£d-l 



/ kJ=0 



(15) 



k.l=0 



£d- 



Hence, this submatrix of ^(M), generated by arbitrary elements with fixed 
Ai , . . . , Ajj^_ 1 , /ii , . . . , /i^M- 1 is ^£^ . This is the last factor of the Kronecker product 
in (13), where c represents one element of the complete previous product. This previous 
product up to the last but one factor is already fulfilling the form of (13) by the induction 
hypothesis. So Lemma 1 implies (13). D 



-1 



mod 



Remark 2. For a slightly loosened version of Lemma 1, i.e. that (h,-, y,) = e 
1, Theorem 1 holds for any pair of biorthogonal bases for ^(M) and ^(M^) and 
even the reverse implication is true: If ,^{M) is of the form (13), then there exist two 
bases fulfilling the slightly loosened version of Lemma 1 and provide an ordering for the 
matrix columns and rows as denoted in (14). The proof is just the reverse steps of the 
proof of Theorem L Hence a pair of biorthogonal bases for ^(M) and '^ {M) fulfilling 
the modified Lemma 1 is necessary and sufficient for the Fourier matrix to fulfill (13). 



4.2 Fast Fourier transform 

Following the approach of Mersereau et. al. [16] or in a more general matter also de- 
scribed by Auslander et al. [2], we can now apply fast Fourier transformation techniques 
by adapting the multivariate Cooley-Tukey-Algorithm. In addition to the latter general 
approach we present a complete algorithm using any biorthogonal bases for the pattern 
^(M) and its dual group ^(M^) from Lemma 1, cf. also Remark 2, and analyse the 
complexity of the algorithm. We also avoid the reindexing mentioned by the former 
authors using the presented basis and their coefficient vectors to arrange the vectors in 
the Fourier transform. 

Given a basis {y j , . . . , y^^^ } of ^ (M) , we can decompose every y G <^ (M) uniquely. 



I.e. 



3!AgEm : y= £Afcy^ 



^(M) 



Using this decomposition every vector b = {by) ^,^r^^ can also be addressed with b 

(^A)AeEM- Let ^ := ^e,.,, ® ^e,.,, ®---® ^e,_, 6 C"X", , 



f- denote the 



Fourier transform with respect to the basis vectors y^, . . . ,yrfj^_i. Using the indexing 
with respect to the position inside the cycles, the Fourier transform reads 



/ ^(0,0,.. .,0) 



= ( 


( 


b 






^M^e. 


^l,2^e. ■ 


■ ^1,.^£, 


= 


^2,l^s, 


^2,2^s, ■ 


■ ^2,.^£, 




^«,l^e. 


^n,2^S. ■ 


■ %i,n'^ed 



\ 



J 



\0.(),...fy-l) 
^...,0,1,0) 



'(0,...,0,l,erf-l) 
V-A2,0) 



\^(ei-Le2-l,...,£rf-l)/ 



(16) 



This enables us to split the computation into two parts: 

First calculating ^e^ ( bn^x^ \ \ for fixed values of Ai , . . . , A^^- 1 • Following 

this blockwise Fourier transform, an interleaved addressing due to fixed e^^ is used 
to compute the multiplications with ^ . This leads to an implementation denoted in 
Algorithm 1. 

Theorem 2 (computational complexity of the FFT on T"^). 

Let M G Z'^^'^, m = I detM| >Qbe given. The Fourier transform b = ^(M)b using (16) 
can be computed with 0{m\ogm) operations. 



as 



Proof. Let any implementation of the one-dimensional FFT be given, i.e. with compu 
tational costs cpvTklogk+0{k) for an input of k coefficients, where e.g. cfft = ^ 
shown in [13]. Then the proof is again by induction over cIm'- The first step is a block- 

wise computation of ^e^ ( ^(Ai,...,a^ ) ) for fixed values of Ai, . . . , A^^-i. These 



«M' 



^%=o 



are n := Ci ■ . . . ■ e^_i Fourier transforms of computation cost CFFT£<i log e^ + 0{£d) each. 
After that, for each fixed < Ajj^ < e^ we have to compute the Fourier transform using 
the Matrix ^, which are e^ transforms, where each of them needs CFFT^logn + 0{n) 
time by induction hypothesis. This leads to 



n (cFFTfirf log £d + 0{£d)) + £d (cFFT^logn + 0{n)) 
CFFT"£rf log £d + 0(m) + CFFT^£d log n + 0{m) 
cfft"^ (log erf + log n) + C>(m) < cfft ''^ log m + 0(m). 



(17) 



D 



10 



1 Fou rierOn Pattern [e, b ] ;= Block[{hatb}, 

2 % e = (e,),-^i: vector containig the elementary divisors of M 

3 % b: vector of input values given as b = {bx)^^T^ 

4 % 

5 % hatb: The Fourier transform b = ^(M)b 

6 If [Length[e] == 1, Return[Fourier[b]]; 

7 % Perform the transforms on blocks of size e^i 

8 Do [ 

9 hatb [[{jUi,...,jUrf^_,}, All]] = Fourier[b[[{jUi,...,jUrfj^_j},AII]]]; 

10 . {Mi.l.ei} {MrfM'l'£rfM-i} 

]; 

12 % Perform transform on the first du-i cycles 

13 % recursively on interleaved blocks 

14 Do [ 

15 % sec denotes c/m^ 1 times the term All 

16 hatb [[ sec, <^]] = FourierOnPattern[{ei, . . . ,e^^ i},hatb[[sec,<^]] 



18 
19 



Return[hatb]; 



Algorithm 1. A fast Fourier transform on ^{M) in Mathematica notation, 
where Fourier denotes any implementation of the one-dimensional Fourier 
transform. 



11 



The proof also reveals the optimisation possibilities of the multivariate Fourier trans- 
form in comparison to the one-dimensional case with the same size m of points: The 
first step consists of n Fourier transforms of the same size e^, where each transform 
acts on a distinct (blockwise) subset of the given data. This means one can use a vec- 
torial SIMD implementation [11] and compute this term in a parallel implementation in 
C)(erflogej). The same holds for the e^ Fourier transforms in the second step, which act 
on distinct (interleaved) subsets on the immediate result, which is in 0{n\ogn), using 
again a vectorial parallel implementation. 

5 Fast wavelet transform 

The same techniques presented in Section 4 can also be applied to the corresponding 
periodic wavelet transform. We first introduce subspaces of L^(T^) that are constructed 
using translates of one function and decompose these into 7" > 2 subspaces. We present 
some properties from [15] concerning these spaces and derive a fast decomposition. We 
also generalise some properties from the one-dimensional case in [21], e.g. the scaling 
property, that is extended by directional information. 

For convenience, we leave out the shift ■ |^ back into the set of congruence class 
representants in the summations onX = ^(M) and X = ^(M^), whenever it is clear 
from context. 

5.1 Translation invariant spaces 

A subspace L C L^(T'^) is called M-invariant if 

VyG^(M):/GL^r(y)/GL. 
Such a subspace can easily be constructed given a function / G L^(T'^). We define 

V^ = span{r(y)/ : y G A(M)} = span{r(y)/ : y G ^(M)} 



as the span of all translates of / (with respect to a matrix M) 
Then for two func 
ily if there exists a 
^(M)a such that 



Then for two functions f,g G L^(T'^') it holds [15, Theorem 3.3] that g G V^ if and 
only if there exists a vector a = («y)y£ j2(m) with its Fourier Transform a = («h)he^(Mr) 



Vh G ^(M^) Vk G Z^ : c^+mtM = «hCh+M^k(/)- (18) 

Then 



12 



If we further look at a decomposition M = JN we see from the definition of r(y), that 
V^ C V^ and it further holds that 

where for a certain order of ^(M^) it follows 

A = (diag (4+jn)he^(N^)),^,,(j,) e C"x- (19) 

Applying this to a set of functions gi,. ..,g\ (jgj ji , whose translates are pairwise orthogo- 
nal, we get a decomposition of V^ into subspaces, cf. [15, Theorem 4.1]. 

5.2 Properties of a multivariate decomposition 

Let J,N G Z'^X'^ be regular matrices and M = JN. We denote the bases of the corre- 
sponding patterns by {xi,...,Xrf^},{yi,...,yrfj^} and {zi,...,Zrfj} for ^{M), ^(N) 
and ^{J). These bases are ordered with respect to the cycle lengths, cf. (8), e.g. y ■ 
corresponds to a cycle of length e^-^+y, i.e. fcy ,■ = mod I <^ ^ = hea-dM+j^ ^ ^ ^• 
We can characterise subspace ,^(N) of ^(M) using the following Lemma. 

vdxd 



Lemma 2 (projection). Let J,N G Z be regular matrices and M = JN. There exists 



a matrix P G Nq'^^ "^ such that for an arbitrary w G ,^(N) we have 



d^ rfjVi 

3Ai G En 3A G Em : w = £ pikyk = Y, ^i^i' (20) 

k=l 1=1 

where A = P/i holds. 

Proof. Equation (20) follows from the fact that ^{N) C ^(M) and decomposing w in 
both bases. Because every basis vector y^ G ^(M) can also be decomposed, we have 

1=1 
This is then applied to 

W = £ HkYk = L ^^ E P/'^'^/ = E E M/t/'/,/tX/ = £ A/X/. 
yt=l /t=l 1=1 1=1 k=l 1=1 

n 



13 



Using this property we can state a generalisation of [2 1 , Theorem 4. 1 .2a] , which deals 
with the translation invariance of subspaces in the one-dimensional case. For any regular 
regular matrix J fulfilling N = J^^M G U^^'^ , the M-invariant space is N-invariant. 
Additionally, we are also able to characterise dimensions and directions of the subspace 
^(M) depending on J as described in the following theorem. In the one-dimensional 
case this would just be the value J G N which is the number of points added per point 
in ,^(N),N G N. Usually one would choose this factor to be 2 to get the classical 
dyadic decomposition scheme of one-dimensional wavelet analysis. In the multivariate 
case however, different matrices J of the same modulus of its determinant exist, which 
enable different extensions of a pattern l^ (N) . These are first stated for some special 
matrices J, i.e. Jj = 1, in the following theorem, which includes all matrices of absolute 
determinant 2. We will discuss the general case after the theorem. 

Theorem 3 (multivariate scaling property). 

Let J,N G Z'^^'^ Z?e regular matrices and M = JN, such that the dimension Jj = 1. 

Denote by £^, £■ , £^, j = 1 , . . . , J, their elementary divisors. Then it holds 

1) forN-^i ^span{y^,...,yj^}, that 

a) Jm = ^n + 1 

b) 3x/G {xi,...,Xrf^} 

N^^zi = Ax/ modi, A G {l,...,ej- 1} andef^ = ej, 

2) for N^^zi G span{y^,. . . ,yrfj^}, that 

a) Jm = dji^ 

b) the decompositions 

N^^zi = J^A/X/, X e Em and N^^zi^Y^l^ky^ M ^ Q"^ 

1=1 k=i 

exist and fulfill 

X - ^Pm, (21) 



where P G Nq'^ ^ ^ is regular 



Proof It holds Jzi =MN-izi gZ^ hence N^^zi G ^(M). If N^^zi ^span{yi, . . . ,y. }, 



then spanjyj , . . . , y^^^ } C spanjy^ , . . . , y^^^ , N ^zi }, which also holds for the spans if we 
restrict the weighted sums to En and En x Ej respectively. 
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From (5) we have for any x G <^(M), the unique decomposition into two elements 

x=y + N-iz|^(j^^, ye^(N)andzG^(J). 
This can be decomposed uniquely using the bases of ^(N) and ^(J), which reads 

k=l 

where the second term is only one summand due to rfj = 1 . Looking at assumption 1 we 
conclude, that the set {yi, . . . ,y^j^,N^^zi} is linear independent, which is a). Statement 
b) follows from the fact that ^ (N) is a subgroup of ^ (M) and hence all cycles of the 
first are also existent in the latter one. This also implies that exactly one basis vector x/ 
spans the new cycle N^^zi . 

Using the same approach, it holds for the second case that N^^zi G spanfy^, . . . ,y^/j^} 
and hence we can decompose 

3!v G En 3!jU G Q^^ 3,9 ^ Ej : N-^zj = £ Vfcy, + £ M^y„ 

k=l k=l 

which is a). Using e^zi G Z"^ and the fact that 

k=l 

we can apply Lemma 2 and see that 

dM ^ dj^ I "^N dM 

3!A gEm:N"^zi = £ A/X/ = ^ £ ^^y^ = ^ £ ^^tJ^P/./tX/, 

/=0 ^d k=0 ^d k=0 1=1 

which is (21). D 

The equality (21) also generalises the notation, that the one-dimensional Fourier 
transform, when the sampling points double in their number (halfening the distance), 
the frequencies double. These one-dimensional Fourier transforms inside the multi- 
variate one can be seen in Theorem 1 and are characterised in Theorem 3 above with 
their directions. In fact even for the first case a cycle, i.e. a one-dimensional Fourier 
transform, gets increased in size, that was e^ = 1 before and hence was not part of the 
basis. 

Theorem 3 can also be generalised to more than one cycle in J by decomposing the 
basis of ^(J) into two distinct parts 81,82 of basis vectors fulfilling the first and the 
second case of the Theorem. Then the dimension gets dM = '^n + l^i I ^nd for each basis 
vector in 82 a scaling property as in (21) holds. 
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5.3 A fast decomposition algorithm 

Given a decomposition M = JN of a regular integral matrix M and the functions 
/,<?!, •••,<?|detj| such that 

|detj| 

^M=©^N' hence ^,= £ ^;,y7^(y)/ e V^ j = 1,.. ., |detj| 

;=1 ye.^(M) 

rf 



we can decompose any 7= Lye^(M) ^yT{y)f G V^ using any pair of bases {hi , . . . , h^^} 
of ^(M^) and {ki, . . . ,krf^} of ^(N^) using the following steps: 

1. Compute a = ^(M)a and bz = ^(M)bz, / = 1,... , |detj| 



2. CalculatePfor the two bases {hi, . . . ,hrf^} of ^(M^) and {ki, . . ..kj^} of ^(N^) 
using Lemma 2. 



3. Calculate for each basis vector of a basis {li , . . . , l^j} of ^(J^) the decomposition 

, (^.,)f=ieEM, j = i,...,Jj 



dm 
N%=£^,,h, 

!=1 



^^(M^) 



and apply these to address N^l G ^(M-^), 1 G ^(J'^) using their coefficient vector 
Xi G Em. 

4. Compute A = P/i Iem for ^^^h h G ^(N^) , which can be reached running through 
all /i G En and compute 

djM = /—-^ £ ^-.h+N^i^+N^i, 7 = l,...,|detJ|, hG^(N^) 

A/|detJ|j^.^(j7-) 

addressing a,b using A + Ailjg and d^-^h using jU. 

5. Perform the inverse Fourier Transform with each d, to get the resulting decompo- 
sition 

|detj| 

r= E E dj,yT{y)gj. 
j=i ye^(N) 

For the special case | det J| = 2, where ^i represents in some sense the low frequent part 
of / and g2 the high frequent part of /, this describes one step in a dyadic fast wavelet 
transform. For any following step in the decomposition, i.e. N = J'O to split V-^^ into 
two spaces, we can apply the steps 2-4 to bi. 
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Given the number m = |detM| of data points as input we obtain the complexity of 
computation for the algorithm. We further note the influence of the dimension d, though 
it is constant with respect to the amount of sampled data points m. The steps 2 and 3 have 
to solve at most | det J| + d linear systems of equations with at most d unknowns each, 
the first step is by (17) 2cFFT"^logm + 0(m), the same without the factor 2 holds for the 
last step, due to |detJ|/7logn < mlogm. Finally, step 4 needs c|detj|m = 0{m) steps, 
where the constant is just depending on the speed of multiplication on that machine. In 
total this is 

3cFFTOTlogm + c|detJ|m + ci|detJ|(i -\-C2{rn-\-d ) 

where c, ci,C2 are constants just depending on speed of multiplication (on a specific 
machine) and they are especially independent of d,m and |detj|. Performing multiple 
decompositions on the same set of data would only affect the second and third term, 
because the Fourier transform would be computed once in the beginning and for every 
wavelet level in total once at the end, which are in total m coefficients. 



6 Example 

6.1 Fourier transform 

As an example for the Fourier transform we look at different patterns ^(M) with the 
same number m = |^(M)| = |detM| of points. The pattern normal form [15, Lemma 
2.5] can be used to look at different classes of matrices that induce the same pattern. 
Restricting the example to the case d = 2, all normal forms have the form 

M=( M , where m = /t/, kj eN+ andO < i <k. (22) 

Furthermore all different cycle lengths Ci , £2, which may occur, are given by all divisors 
/ of k. 

Given any implementation of the one-dimensional Fourier transform, the Algorithm 1 
can be easily implemented. For any of the recursive calls in line 16, the set of data a 
subroutine is working on, is distinct from all the other. Hence, the calls can also be 
parallelized. 

We compare the implementation of a serial and a parallel algorithm of the multivariate 
Fourier transform, where the first one is compared to the usual one-dimensional Fourier 
transform on the same amount of data, i.e. an m-dimensional vector. The algorithms 
are all using the implementation of Fourier given by Mathematica 8 running on an Intel 
Pentium 4 Core 2 Quad with 2 Ghz each and 4 GB memory. 

Choosing k = I = 2048 and a randomly generated set of m = iP data values, the 
computational times for all divisors i of k are depicted in Table 1 . The factor denoted 
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i 


cycles 


serial (sec.) 


factor 


parallel (sec.) 


gain 


1 


(4194304) 


0.266664 


1.02242 


0.270647 


0.98528 


2 


(2 2097152) 


0.471722 


1.80864 


0.343969 


1.37141 


4 


(4 1048576) 


0.468220 


1.79522 


0.299426 


1.56373 


8 


(8 524288) 


0.457126 


1.75268 


0.282309 


1.61924 


16 


(16 262144) 


0.457321 


1.75342 


0.286432 


1.59661 


32 


(32 131072) 


0.469080 


1.79851 


0.299551 


1.56595 


64 


(64 65536) 


0.466902 


1.79016 


0.315000 


1.48223 


128 


fl28 32768) 


0.572192 


2.19386 


0.412570 


1.38690 


256 


(256 16384) 


0.920226 


3.52826 


0.741036 


1.24181 


512 


(512 8192) 


1.144126 


4.38672 


0.530716 


2.15582 


1024 


(1024 4096) 


0.949403 


3.64013 


0.424802 


2.23493 





(2048 2048) 


0.907286 


3.47865 


0.411705 


2.20373 



Table 1. For all possible cycles of M from (22), where k = l = 2^^ we compare 
a serial and a parallel implementation of the two-dimensional FFT. The first 
is compared to the one-dimensional FFT of m = | detM| = 2^^ points of data 
using Fourier in Mathematica. This took 0.260816 seconds, which leads to the 
factor (col. 4). The table lists all possible elementary divisors / of k, which 
corresponds to the obtained cycles. The gain represents the factor, the parallel 
implementation gains — on 4 cores — compared to the serial one. All times are 
obtained taking the mean value of 50 measurements using AbsoluteTiming. 
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in the column after the serial one is obtained by comparison to the one-dimensional 
Fourier transform. The multidimensional data is obtained by sampling along the cycles 
of ^(M), which represent the sampling directions. This data is then used to perform 
the lattice Fourier transform and measuring its serial and parallel computation time. The 
last column denotes the parallel gain. 

While the serial implementation is equal to the one-dimensional case for / = 1, it is 
slower for the cases, where the Fourier transform has to be applied to rows and columns. 
The parallel implementation is able to gain a factor of 1.2 to 2.8 on 4 cores. 

The lattice Fourier transform reduces with a change of basis to the usual multivari- 
ate Fourier transform. Further investigations of the parallelization of the multivariate 
Fourier transform can be found e.g. based on the FFTW in [19]. 

6.2 Wavelet transform 

To extend the example from the previous subsection, we look at the Dirichlet type 
wavelets for different decompositions M = JN, where 

Let ^(M) = A(M) n [-i, ^y and ^(M) = M^(M) be full sets of congruence class 
representants of A(M) and Z'^. We denote by r(k) = rM(k),k G Z^ the number of 
superficial hyperplanes of the parallelepiped M^[— ^, jY apoint h = k|f^(-jy[r) G ^(M^) 
is lying on, i.e. 



r(k)=#|;:|h^M-i|, = i, h=k^Mr)|. 



The Dirichlet kernel 9m : T'' — )■ M is given by its Fourier coefficients 

r(k) 
2 

else. 



Ck(<?>M) = <^ ^ ^ 2'2J ' kGZ^. 



The translates T{y)(pM, y G ^(M) are linear independent and orthonormal [15, Theo- 
rem 6.4] and it holds ^ G V^ due to 



Ck(9N) =^hCk(9M), kGZ'', h=ky(j^T), where 



else. 
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Figure 2. The centered box splines 5s (left) and 5>j/ (right) used as test func- 
tions for the Dirichlet type wavelet decomposition. 



see (18) and [15, Lemma 6.5]. The corresponding wavelet i/^N fulfilling Vj 
V^^ is given by 



<Pm 

M 



C' 



Ck(V^N) 







g|<^(M7') 



,-2mk^N"'y 



ifM-^kG 
else, 



'2'2J ' 



where y G ,^(J)\{0} and g G ^(J^)\{0} are uniquely determined due to |detj| 
for each J G ^ . 



As test functions we choose two centered box splines 5s,5»p 
?\ _ /;r f f 



— ;■ 



;r 



^ 



;r 



with S = 
, see [3] and Fig. 2. The function Bz is a piece- 



wise linear polynomial, 5ij/ a piecewise cubic polynomial. We choose M : 



512 






512 



and the function fv G V^ is obtained by sampling B^ at the points 2;ry, y G ^(M). 
From these samples we perform a change of basis from the interpolatory or Lagrangian 
basis into the basis of orthonormal translates {7'(y)<?>M}ye^(M) by implementing the 
well known change of basis, which is for this case stated e.g. in [15, Corollary 3.7]. 

We perform the first step of the wavelet transform, i.e. we decompose f = fv +fw, 
where fy G V^ and fw G V^^ . The function /ly is in the 17th wavelet space of a dyadic 
pyramid of wavelet spaces due to | detN| =2^^. 

Depending on the matrix J we obtain different directional information in the wavelet 
space for the box spline Bz'. The lines of discontinuity of the first derivative are divided 
into three sets: The discontinuities along the diagonal can be obtained using J^ (see 
Fig. 3(c)), while discontinuities parallel to the axes are depicted in the wavelet spaces 
generated by i^ and J^, (see Figs. 3(a) and 3(b)). The functions /ly are shown in their 
absolute value to emphasise the nonzero parts of these wavelet spaces. The oscillations 
seen as grey values between the lines or to the rim are due to the Dirichlet type wavelets. 
The same holds for the diagonal discontinuities that are also visible in the first two 
decompositions. 
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A 



A 



/ 



V 



V 



/ 



(a) l/wl using J;t 



(b) l/wl using Jy 



(c) |/h/| using Jrf 



Figure 3. The three possible functions \fw\ as wavelet spaces of 5^ drawn as 
a contour plot, where white indicates zero and black the highest value. 



// 



/j 



V 



(a) |/w| using J;c 



(b) |/w| using Jy 



(c) I/H'lusingjrf 



Figure 4. The three possible functions |/iy | as wavelet spaces of 5>j< drawn as 
a contour plot, where white indicates zero and black the highest value. 
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For the second example we apply the Dirichlet type wavelets to B^i (see Fig. 4). 
The wavelet spaces contain the discontinuities of the third derivative. Here, we see the 
same oscillations as in the previous case, but the decomposition with respect to certain 
directions is more clear. The union of all three wavelet spaces represents the set of lines, 
on which the third derivative is discontinuous along the directional derivative orthogonal 
to the corresponding hne. 
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